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Abstract 



First principles calculations based on density functional theory are having 
an increasing impact on our understanding of molecule-surface interactions. 
For example, calculations of the multi-dimensional potential energy surface 
have provided considerable insight into the dynamics of dissociation processes. 
However, these calculations using a plane- wave basis set are very compute 
expensive if they are to be fully converged with respect to the plane-wave 
energy cutoff, k-point sampling, supercell size, slab thickness, etc. Because 
of this, in this study, we have implemented a mixed-basis approach which 
uses pseudo-atomic orbitals and a few low-energy plane waves as the basis set 
within a density functional, pseudopotential calculation. We show that the 
method offers a computationally cheap but accurate alternative. The energy 
barrier for hydrogen dissociation on Cu(lll) is calculated as an example. 

Keywords: Chemisorption, Copper, Density functional calculations, Ab-initio 
quantum chemical methods and calculations, Hydrogen, Low index single 
crystal surfaces, Models of surface chemical reactions 
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I. INTRODUCTION AND COMPUTATIONAL METHOD 



In recent years, the application of first principles electronic structure methods to surface 
systems has increased significantly as a result of improvements in algorithms and enhanced 
computational speed. The system size and complexity which can be analysed has increased 
by an order of magnitude following the pioneering work of Car and Parrinello M, who 
introduced an iterative minimisation of the total energy based on wavefunction improvements 
at each iteration. The construction of the initial wavefunction is clearly important for 
efficiency in such an approach. 

Among several different approaches, there are two simple and natural choices of basis 
set for the expansion of electron wavefunctions: atomic orbitals and plane waves. On the 
negative side, atomic orbital methods have difficulties in representing the wavefunctions 
and potential in interstitial and vacuum regions while plane wave expansions are expensive 
for representing localised atomic character, for example 3d wavefunctions. Nevertheless, 
plane-wave basis sets are in most common use since they are simple, independent of atomic 
positions, fast Fourier transformation (FFT) methods can be applied readily, and accuracy 
can be systematically improved by including additional plane waves with higher energy cut- 
offs. Although atomic orbitals are more physical it is difficult to represent uniform charge 
density, as in the vacuum region of a surface, with atom-centered, localised orbitals. On 
the other hand, plane-wave basis sets are also inefficient in a surface calculation using a 
slab geometry, since one needs as many plane waves for the vacuum as for the solid region. 
Therefore, a combination of the important properties of plane waves with atomic orbitals in 
a mixed basis may give a convenient and efficient representation, especially for systems which 
include both highly localised (atomic-like) and delocalised (plane-wave-like) components. 

In this study, we have implemented a mixed-basis approach || which uses pseudo-atomic 
orbitals and a few low-energy plane waves as the basis set within a density functional, 
pseudopotential calculation. A similar approach has been described by Neugebauer and 
Van de Walle ||, but they focussed on providing a better starting wavefunction for a full 
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plane-wave calculation. Here, we investigate whether the mixed-basis calculation might be 
used in its own right for calculations of the potential energy surface for molecule-surface 
systems. 

The Kohn-Sham eigenf unctions are expanded as 



where a is the band index, \x is a combined index which labels the orbitals and atomic sites, 
a and b are coefficients of the pseudo-atomic orbitals and plane waves, respectively, and Q 
is the volume of the unit cell. x», is the Bloch sum formed from pseudo-atomic orbitals as 



where m labels the orbitals, the Ri are the lattice vectors, the t\ are the atomic coordinates, 
and m are pseudo- atomic orbitals. In practice, we use a plane- wave expansion for x*i(^)> 
and exactly the same FFT grid as in a full plane-wave calculation. There are therefore two 
plane-wave energy cut-offs to be considered in the mixed-basis calculation. The larger one 
is the cut-off used in the representation of x^if) and is the same as would be used in a 
full plane-wave calculation. The smaller one is the cut-off for the extra, low-energy plane 
waves which appear in the second term of Eq. |l|. This plane- wave representation of Xnif) 
makes the calculation of the charge density, the kinetic energy, multicentre integrals, and 
the contribution from non-local pseudopotentials straightforward. 

Solving the Schrodinger equation then reduces to solving the secular equation 



The overlap matrix elements are given by (with reference to the partition of ip in Eq.p 
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where v = (n, j) and I^{k) is the Fourier integral of the pseudo-atomic orbital, 

If{k) = j=J dre-^^ m {f). (7) 



Similarly, the Hamiltonian matrix elements are based on a plane-wave representation 



H S3 , = J= / dre-^+^He 1 ^ 3 ^ (8) 



1 

fc + G| 2 «5 gd > + V^(G - G) + V NL ({k + G), {k + G')). (9) 



Then 



H»g = Y,^ I ?$)H, (} (10) 
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The local part of the potential, Vi oca i, in Eq. |9] contains the Hartree and exchange-correlation 
potentials, as well as the local part of the pseudopotential. In practice, only the Hartree 
and exchange-correlation contributions need be re-calculated through the self-consistency 
cycle — the pseudopotential (both local and non-local, Vnl, parts) and kinetic energy matrix 
elements are calculated only at the first iteration. Self-consistency is achieved by a combi- 
nation of Kerker charge density mixing and a modified Broyden method @|||. The initial 
charge density is constructed from overlapping, atomic pseudo-charge densities. 

Diagonalization of Eq. || is acceptable, since there are at most 9 orbitals (s,p,d) for each 
atom, and typically 10 to 20 additional plane waves per atom (see below). This results in 
a matrix size less than 10 3 x 10 3 , even for a moderately large system, compared to between 
10 4 to 10 5 for a pure plane-wave expansion. In fact, for the H 2 /Cu system considered below, 
the most expensive part of the whole calculation is the construction of the Hamiltonian 
matrix. Tests on this system show that the mixed-basis method is 6 to 8 times faster per 
iteration than our pure plane wave code (as described in 0). In addition, it typically requires 
fewer than half as many iterations to converge, and so provides a significant improvement 
in computational speed. A full analysis of the timing and scaling of the computation with 
respect to the important calculational parameters will be presented elsewhere. 



II. CALCULATIONS AND RESULTS 



We have carried out a careful benchmarking analysis of the mixed-basis approach by 
increasing the number of additional low-energy plane-waves and comparing the results with 
an exactly equivalent calculation based on plane waves only. The system chosen for this 
comparison is hydrogen dissociation on Cu(lll), which is a model system for the experi- 
mental (eg [0]) and theoretical (eg PHnj) study of dissociation dynamics. In particular, 
we concentrate on the value of the minimum energy barrier, which is known to occur for 
the geometry in which the molecular axis is parallel to the surface plane and with the H 2 
molecule over a bridge site with the H atoms pointing towards neighbouring hollow sites |§ . 
In all calculations the transition state is taken to be where the H2 molecule is 1.2A above 
the top- layer Cu atoms, with a bond length of 1.1 A |10 . 



As a preliminary study, we have examined some properties of bulk fee Cu and the H 2 

pseu- 



molecule within the mixed-basis scheme. A semi-relativistic, Troullier-Martins [11 



dopotential (with associated pseudo-atomic orbitals) is used to describe the Cu atoms, and 
hydrogen is described by the full Coulombic potential, with localised Is and 2p orbitals. 
For Cu, the irreducible wedge of the fee Brillouin zone is sampled with 28 Appoints and the 
H 2 molecule is calculated in the same cell as in the H 2 /Cu(lll) system described below. 
The calculated lattice parameter and bulk modulus for Cu; and bond length and vibrational 
frequency (estimated from a harmonic fit about the equilibrium bond length) for H 2 are pre- 
sented in Tables I and II respectively, as a function of the cut-off energy for the low-energy 
plane waves. It can be seen that for these two very different systems the mixed-basis method 
converges rapidly to the full plane- wave result. 

Technical details of the full H 2 /Cu(lll) calculation are as follows. The substrate is 
modelled by a 5 layer, rigid Cu slab with the experimental lattice parameter of 3.6lA. There 
are 3 Cu atoms per layer within a x y/3 geometry and 4 layers of vacuum separate 
the slabs. An H 2 molecule is placed on only one side of the slab. The surface Brillouin 
zone is sampled by 54 special /c-points (15 in the irreducible wedge), and the Fermi surface 
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is broadened with a 0.25eV smearing function, with the total energy being extrapolated 
to zero temperature. A cut-off of 800eV is used for the pure plane-wave calculation, and 
for the representing the localised orbitals in the mixed-basis method. These calculational 
parameters should provide well-converged results for the energy barrier ||10| . The local 
density approximation (LDA) is used in the self-consistency cycle of the calculations, and 
the generalised gradient approximation (GGA) (see eg |12|) is incorporated by calculating 
the exchange-correlation energy of the LDA density. 

Six values of the minimum barrier have been calculated. Five use the mixed-basis ap- 
proach with different numbers of additional plane waves, corresponding to energy cut-offs 
of OeV (ie pure pseudo-atomic orbitals), 20eV, 40eV, 60eV and 80eV. The typical number 
of additional plane waves are 0, 65, 180, 330 and 520 respectively for these cut-offs. The 
sixth value corresponds to the full plane-wave calculation. Computed values of the minimum 
energy barrier are given in Table fTJ. It can be seen that the mixed-basis method provides 
an accurate value for the barrier height (within 10 meV) once plane waves with energies up 
to 40 eV are included. 



III. CONCLUSIONS 

The mixed-basis method we have presented appears to provide an accurate and com- 
putationally cheap alternative to full plane-wave methods for first principles calculations of 
surface systems. Although we have discussed only a single case in detail, tests on a variety 
of systems have confirmed the generality of this approach, with structural parameters being 
accurate to of order O.OlA and energy differences accurate to a few tens of meV. Only total 
energies have been discussed here, but it is also straightforward to calculate forces within 
the mixed-basis scheme using the general force expressions given in ||13| . This holds out 
the prospect of being able to perform, with high computational efficiency, a full structural 
relaxation of a system using the mixed-basis approach. If required, the calculation can then 
be "finished off" using a full plane-wave expansion, in which case the mixed-basis method 
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will provide an accurate starting structure and highly optimised initial wavef unctions. 
This work is supported by the UK Engineering and Physical Sciences Research Council. 
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TABLES 



Energy cut-off (eV) Lattice constant (A) Bulk Modulus (GPa) 

AO 3.59 165 

20 3.57 180 

40 3.57 177 

60 3.58 173 

80 3.58 175 

PW 3.57 178 



TABLE I. Calculated LDA lattice constant and bulk modulus for fee Cu, as a function of the 
cut-off energy for additional plane waves in the mixed-basis approach. AO represents the pure 
pseudo-atomic orbital limit (ie eV cut-off) and PW the full plane-wave limit (with plane waves 
up to 800 eV). The Murnaghan equation of state is used. 



Energy cut-off (eV) Bond length (A) Vibrational quantum (meV) 

AO 0.837 535 

20 0.781 489 

40 0.761 471 

60 0.756 467 

80 0.754 464 

PW 0.753 462 



TABLE II. Calculated GGA bond length and vibrational quantum (hto) for H2. Symbols are 
as in Table I. 
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Energy cut-off (eV) 



Barrier Height (eV) 



AO 
20 
40 
60 
80 

PW 



-0.854 
0.520 
0.606 
0.614 
0.620 
0.612 



TABLE III. GGA minimum energy barrier for H2 dissociation on Cu(lll). Symbols are as in 
Table I. 
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